! Definition of
! 
! NodeBoundaryMeshLabel(:), nodeInterface
! sendCountsVPoints(:), recvCountsVPoints(:)
! rDisplsVPoints(:), sDisplsVPoints(:)
!
! sendVPoints9(:), recvVPoints9(:)            
! sDisplsVPoints9(:), rDisplsVPoints9(:)
!
! allocation of
! sendVPoints(nodeInterface),recvVPoints(nodeInterface)
! sendCountsVPoints9(:), recvCountsVPoints9(:)
subroutine CalcNodeInterface()
use Data_cfoam
use Data_Mesh, only: nPoints
implicit none
include 'mpif.h'

    integer,allocatable :: CountsPoints(:),DisplsPoints(:)
    integer,allocatable :: LocalOrder(:)
    integer,allocatable :: sendpoint(:),recpoint(:)
    integer :: iProc, iPoint
    integer :: i, j, i1,i2, j1,j2
    integer :: statusp(mpi_status_size),ifail
    character :: sfile(3)
    
    nodeInterface = 0
    sDisplsVPoints = 0
    rDisplsVPoints = 0
    recvCountsVPoints = 0
    sendCountsVPoints = 0
    sDisplsVPoints9 = 0
    rDisplsVPoints9 = 0
    recvCountsVPoints9 = 0
    sendCountsVPoints9 = 0
    
    if(nump==1)return
    
    allocate(CountsPoints(nump),DisplsPoints(nump+1))
    CountsPoints = 0
    
    call MPI_AllGather(nPoints,1,MPI_INTEGER,CountsPoints,1,MPI_INTEGER,MPI_COMM_WORLD,ierr)
    iPoint=maxval(CountsPoints)

    DisplsPoints(1) = 0
    do iProc = 1,nump
      DisplsPoints(iProc+1) = DisplsPoints(iProc) + CountsPoints(iProc)
    end do
    
!    do iProc = 1,nump
!        write(2,*) '@@@ np=', rank, ' nCountsP(', iProc, ')=', CountsPoints(iProc), &
!            ' DisplsP=', DisplsPoints(iProc+1)
!    end do
!    write(2,*) '------------------------------'
    
    allocate(LocalOrder(nPoints))
    allocate(sendpoint(iPoint),recpoint(iPoint))

    do i = 1,nPoints
      LocalOrder(i) = DisplsPoints(rank+1)+i
    end do

    call Reodering(pointProcAddressing,LocalOrder,0,nPoints,nPoints,ifail)

    if(ifail/=0)STOP '@FATAL ERROR in Reordering!'

! put sendrec stuff here
    do i = 1,nPoints
      LocalOrder(i) = LocalOrder(i)-DisplsPoints(rank+1)
      sendpoint(i) = pointProcAddressing(LocalOrder(i))
    end do

! determine value of nodeInterface to enable dynamic memory allocation of NodeBoundaryMeshLabel
    do iProc = 0,nump-1
    IF (iProc/=rank)THEN
    call MPI_SENDRECV(sendpoint(1), nPoints, MPI_INTEGER, iProc, 0, &
        recpoint(1), CountsPoints(iProc+1), MPI_INTEGER, iProc, 0, &
        MPI_COMM_WORLD,statusp, ierr)

    do i = 1,nPoints
    do j = 1,CountsPoints(iProc+1)
        IF(sendpoint(i)==recpoint(j))THEN
        nodeInterface = nodeInterface + 1
        ENDIF
    end do
    end do
    ENDIF
    end do

    allocate(NodeBoundaryMeshLabel(nodeInterface))

! reset nodeInterface to enable buildup of NodeBoundaryMeshLabel
    nodeInterface = 0

    do iProc = 0,nump-1
    IF (iProc/=rank)THEN
    call MPI_SENDRECV(sendpoint(1), nPoints, MPI_INTEGER, iProc, 0, &
        recpoint(1), CountsPoints(iProc+1), MPI_INTEGER, iProc, 0, &
        MPI_COMM_WORLD,statusp, ierr)
!        write(*,*) rank,'sendrec from', iProc

    do i = 1,nPoints
    do j = 1,CountsPoints(iProc+1)
        IF(sendpoint(i)==recpoint(j))THEN
        nodeInterface = nodeInterface + 1
        NodeBoundaryMeshLabel(nodeInterface)=LocalOrder(i)
        ENDIF
    end do
    end do
    ENDIF

    sDisplsVPoints(iProc+2:nump) = nodeInterface
    recvCountsVPoints(iProc+1) = nodeInterface - sDisplsVPoints(iProc+1)
    if(recvCountsVPoints(iProc+1)/=0)write(*,*) rank,'neigh',iProc
    end do

    do i = 1,nump
      sendCountsVPoints(i) = recvCountsVPoints(i)
      if(recvCountsVPoints(i).eq.0) sDisplsVPoints(i) = 0
      rDisplsVPoints(i) = sDisplsVPoints(i)  
      sDisplsVPoints9(i) = 9 * sDisplsVPoints(i)
      rDisplsVPoints9(i) = sDisplsVPoints9(i)
      recvCountsVPoints9(i) = 9 * recvCountsVPoints(i)
      sendCountsVPoints9(i) = recvCountsVPoints9(i)
!	    write(2,*) '@@@@ rank=',rank, ' recvCounts=',recvCountsVPoints(i),' sDispls=',sDisplsVPoints(i)
    end do
	
!    write(2,*) '------------------------------'
    
    allocate(sendVPoints(nodeInterface),recvVPoints(nodeInterface))
    allocate(sendVPoints9(9*nodeInterface),recvVPoints9(9*nodeInterface))

!    do i = 1,nodeInterface
!        write(2,'(a5,i3,a3,i3,a6,i5,a7,3f10.5)') 'rank=', rank, ' i=', i, ' Node=', NodeBoundaryMeshLabel(i), ' Coord=', meshPoints(1:3,NodeBoundaryMeshLabel(i))
!    end do

!    write(*,*) '@@@ nodeInterface=', nodeInterface

    deallocate(CountsPoints,DisplsPoints)
    deallocate(LocalOrder)
    deallocate(sendpoint,recpoint)

end subroutine CalcNodeInterface
! definition of transposition from min to max, from i1+1 to i2
! in: iArray - array 1:n
!     i1,i2 - interval in iArray for transposition, i1<=i2<=n
! out: Order - transposition: iArray(Order(i))<iArray(Order(j)), if i<j
!      ierr = 0 if Reodering have done successfully
subroutine Reodering(iArray,Order,i1,i2,n, ierr)
implicit none
    integer,intent(in) :: iArray(1:n),i1,i2,n
    integer,intent(inout) :: Order(1:n)
    integer,intent(out) :: ierr

    integer :: i,j, length, aux
    integer :: iLocalArray(i2-i1),iLocalOrder(i2-i1)

    ierr = 0
    if(i1.gt.i2.or.i2.gt.n)then
        ierr = 1
        return
    end if

    length = i2-i1
    iLocalArray(1:length) = iArray(i1+1:i2)
    iLocalOrder(1:length) = Order(i1+1:i2)

    do i = 1,length-1
        do j = i+1,length
            if(iLocalArray(i).gt.iLocalArray(j))then
                aux = iLocalArray(i)
                iLocalArray(i) = iLocalArray(j)
                iLocalArray(j) = aux

                aux = iLocalOrder(i)
                iLocalOrder(i) = iLocalOrder(j)
                iLocalOrder(j) = aux
            end if
        end do
    end do

    Order(i1+1:i2) = iLocalOrder(1:length)
end subroutine Reodering


